# Replication file
# Election models in "Fixed Effects and Post-Treatment Bias in Legacy Studies"

rm(list=ls())
set.seed(12435)
library(stargazer)

load("elections_fes.RData")

##
## Analysis
##
elect_data_full <- elect_data
elect_data70 <- subset(elect_data, distance2<7)


####
#### Approach 1: contemporary states in first stage only
####
# a) AFD, full sample
mod01_a <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(NAME_1), elect_data)
mod01_b <- lm(I(AfDshare - coef(mod01_a)["foreign_rate"]*foreign_rate 
                - coef(mod01_a)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden, elect_data)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data[sample(1:nrow(elect_data), replace=TRUE),]
  boot1 <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(NAME_1), databoot)
  boot2 <- lm(I(AfDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden, databoot)
  gboot[b,] <- coef(boot2)
  if (b%%250==0) print(b)
}
SEs1 <- apply(gboot,2,sd)
names(SEs1) <- names(coef(boot2))

# b) AFD, 70km
mod01_c <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(NAME_1), elect_data70)
mod01_d <- lm(I(AfDshare - coef(mod01_c)["foreign_rate"]*foreign_rate 
                - coef(mod01_c)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden, elect_data70)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data70[sample(1:nrow(elect_data70), replace=TRUE),]
  boot1 <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(NAME_1), databoot)
  boot2 <- lm(I(AfDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden, databoot)
  gboot[b,] <- coef(boot2)
  if (b%%250==0) print(b)
}
SEs2 <- apply(gboot,2,sd)
names(SEs2) <- names(coef(boot2))

# c) AFD+NPD, full sample
mod01_e <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(NAME_1), elect_data)
mod01_f <- lm(I(AfDNPDshare - coef(mod01_e)["foreign_rate"]*foreign_rate 
                - coef(mod01_e)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden, elect_data)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data[sample(1:nrow(elect_data), replace=TRUE),]
  boot1 <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(NAME_1), databoot)
  boot2 <- lm(I(AfDNPDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden, databoot)
  gboot[b,] <- coef(boot2)
  if (b%%250==0) print(b)
}
SEs3 <- apply(gboot,2,sd)
names(SEs3) <- names(coef(boot2))

# d) AFD+NPD, 70km
mod01_g <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(NAME_1), elect_data70)
mod01_h <- lm(I(AfDNPDshare - coef(mod01_g)["foreign_rate"]*foreign_rate 
                - coef(mod01_g)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden, elect_data70)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data70[sample(1:nrow(elect_data70), replace=TRUE),]
  boot1 <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(NAME_1), databoot)
  boot2 <- lm(I(AfDNPDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden, databoot)
  gboot[b,] <- coef(boot2)
  if (b%%250==0) print(b)
}
SEs4 <- apply(gboot,2,sd)
names(SEs4) <- names(coef(boot2))


###
### TABLE SM3.1
###
stargazer(mod01_b, mod01_d, mod01_f, mod01_h,
          se=list(SEs1,SEs2,SEs3,SEs4),
          star.cutoffs = c(0.05, 0.01, NA), 
          column.sep.width="1pt", digits=3, no.space=T,
          dep.var.labels.include = FALSE)

cfs_currentstates_1st <- c(coef(mod01_b)[2], coef(mod01_d)[2], coef(mod01_f)[2], coef(mod01_h)[2])
ses_currentstates_1st <- c(SEs1[2],SEs2[2],SEs3[2],SEs4[2])



####
#### Approach 2: Weimar states based on centroids
####
# a) AFD, full sample
mod01_a <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid), elect_data)
mod01_b <- lm(I(AfDshare - coef(mod01_a)["foreign_rate"]*foreign_rate 
                - coef(mod01_a)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), elect_data)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data[sample(1:nrow(elect_data), replace=TRUE),]
  boot1 <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid), databoot)
  boot2 <- lm(I(AfDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs1 <- apply(gboot,2,sd)

# b) AFD, 70km
mod01_c <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid), elect_data70)
mod01_d <- lm(I(AfDshare - coef(mod01_c)["foreign_rate"]*foreign_rate 
                - coef(mod01_c)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), elect_data70)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data70[sample(1:nrow(elect_data70), replace=TRUE),]
  boot1 <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid), databoot)
  boot2 <- lm(I(AfDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs2 <- apply(gboot,2,sd)

# c) AFD+NPD, full sample
mod01_e <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid), elect_data)
mod01_f <- lm(I(AfDNPDshare - coef(mod01_e)["foreign_rate"]*foreign_rate 
                - coef(mod01_e)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), elect_data)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data[sample(1:nrow(elect_data), replace=TRUE),]
  boot1 <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid), databoot)
  boot2 <- lm(I(AfDNPDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs3 <- apply(gboot,2,sd)

# d) AFD+NPD, 70km
mod01_g <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid), elect_data70)
mod01_h <- lm(I(AfDNPDshare - coef(mod01_g)["foreign_rate"]*foreign_rate 
                - coef(mod01_g)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), elect_data70)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data70[sample(1:nrow(elect_data70), replace=TRUE),]
  boot1 <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid), databoot)
  boot2 <- lm(I(AfDNPDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs4 <- apply(gboot,2,sd)

###
### TABLE SM3.2
###
stargazer(mod01_b, mod01_d, mod01_f, mod01_h,
          se=list(SEs1,SEs2,SEs3,SEs4),
          star.cutoffs = c(0.05, 0.01, NA), 
          column.sep.width="1pt", digits=3, no.space=T,
          dep.var.labels.include = FALSE)

cfs_weimar_centroid <- c(coef(mod01_b)[2], coef(mod01_d)[2], coef(mod01_f)[2], coef(mod01_h)[2])
ses_weimar_centroid <- c(SEs1[2],SEs2[2],SEs3[2],SEs4[2])



####
#### Approach 3: contemporary states in first stage AND Weimar states based on centroids
####
# a) AFD, full sample
mod01_a <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid) + as.factor(NAME_1), elect_data)
mod01_b <- lm(I(AfDshare - coef(mod01_a)["foreign_rate"]*foreign_rate 
                - coef(mod01_a)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), elect_data)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data[sample(1:nrow(elect_data), replace=TRUE),]
  boot1 <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid) + as.factor(NAME_1), databoot)
  boot2 <- lm(I(AfDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs1 <- apply(gboot,2,sd)

# b) AFD, 70km
mod01_c <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid) + as.factor(NAME_1), elect_data70)
mod01_d <- lm(I(AfDshare - coef(mod01_c)["foreign_rate"]*foreign_rate 
                - coef(mod01_c)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), elect_data70)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data70[sample(1:nrow(elect_data70), replace=TRUE),]
  boot1 <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid) + as.factor(NAME_1), databoot)
  boot2 <- lm(I(AfDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs2 <- apply(gboot,2,sd)

# c) AFD+NPD, full sample
mod01_e <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid) + as.factor(NAME_1), elect_data)
mod01_f <- lm(I(AfDNPDshare - coef(mod01_e)["foreign_rate"]*foreign_rate 
                - coef(mod01_e)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), elect_data)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data[sample(1:nrow(elect_data), replace=TRUE),]
  boot1 <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid) + as.factor(NAME_1), databoot)
  boot2 <- lm(I(AfDNPDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs3 <- apply(gboot,2,sd)

# d) AFD+NPD, 70km
mod01_g <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid) + as.factor(NAME_1), elect_data70)
mod01_h <- lm(I(AfDNPDshare - coef(mod01_g)["foreign_rate"]*foreign_rate 
                - coef(mod01_g)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), elect_data70)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data70[sample(1:nrow(elect_data70), replace=TRUE),]
  boot1 <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid) + as.factor(NAME_1), databoot)
  boot2 <- lm(I(AfDNPDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs4 <- apply(gboot,2,sd)

###
### TABLE SM3.3
###
stargazer(mod01_b, mod01_d, mod01_f, mod01_h,
          se=list(SEs1,SEs2,SEs3,SEs4),
          star.cutoffs = c(0.05, 0.01, NA), 
          column.sep.width="1pt", digits=3, no.space=T,
          dep.var.labels.include = FALSE)

cfs_weimar_contemp <- c(coef(mod01_b)[2], coef(mod01_d)[2], coef(mod01_f)[2], coef(mod01_h)[2])
ses_weimar_contemp <- c(SEs1[2],SEs2[2],SEs3[2],SEs4[2])



####
#### Plots for FIGURE 4
####

#### AFD full sample
coeffs <- c(-.081, cfs_currentstates_1st[1], cfs_weimar_centroid[1], cfs_weimar_contemp[1])
sts <- c(.016, ses_currentstates_1st[1], ses_weimar_centroid[1], ses_weimar_contemp[1])
pdf('elections_plot1.pdf', width=6.25, height=4.2)
par(mar=c(2.5,8,1,2))
plot(seq(-.15,.05, length.out = 2), c(.5,4.5), type='n',yaxt='n', ylab='',
     xlab='',cex.lab=.9)
abline(v=0,lty=2)
pos <- 1:4
axis(side= 2,at=rev(pos),labels=c('Baseline', 'Current state FEs\nin 1st stage',
                                  'Weimar state FEs\n(centroid interpolation)',
                                  'Current state FEs\n+ Weimar state FEs'),las=2,tick=F, cex.axis=.8)
points(coeffs,rev(pos),pch=20)
for (i in 1:4){
  segments(coeffs[i]-1.96*sts[i],rev(pos)[i],
           coeffs[i]+1.96*sts[i])
  segments(coeffs[i]-1.65*sts[i],rev(pos)[i],
           coeffs[i]+1.65*sts[i], lwd=2)
}
dev.off()

#### AFD 70km
coeffs <- c(-.159, cfs_currentstates_1st[2], cfs_weimar_centroid[2], cfs_weimar_contemp[2])
sts <- c(.062, ses_currentstates_1st[2], ses_weimar_centroid[2], ses_weimar_contemp[2])
pdf('elections_plot2.pdf', width=4.3, height=4.2)
par(mar=c(2.5,2,1,2))
plot(seq(-.35,.1, length.out = 2), c(.5,4.5), type='n',yaxt='n', ylab='',
     xlab='',cex.lab=.9)
abline(v=0,lty=2)
pos <- 1:4
points(coeffs,rev(pos),pch=20)
for (i in 1:4){
  segments(coeffs[i]-1.96*sts[i],rev(pos)[i],
           coeffs[i]+1.96*sts[i])
  segments(coeffs[i]-1.65*sts[i],rev(pos)[i],
           coeffs[i]+1.65*sts[i], lwd=2)
}
dev.off()

#### AFD + NPD full sample
coeffs <- c(-.092, cfs_currentstates_1st[3], cfs_weimar_centroid[3], cfs_weimar_contemp[3])
sts <- c(.017, ses_currentstates_1st[3], ses_weimar_centroid[3], ses_weimar_contemp[3])
pdf('elections_plot3.pdf', width=6.25, height=4.2)
par(mar=c(2.5,8,1,2))
plot(seq(-.15,.05, length.out = 2), c(.5,4.5), type='n',yaxt='n', ylab='',
     xlab='',cex.lab=.9)
abline(v=0,lty=2)
pos <- 1:4
axis(side= 2,at=rev(pos),labels=c('Baseline', 'Current state FEs\nin 1st stage',
                                  'Weimar state FEs\n(centroid interpolation)',
                                  'Current state FEs\n+ Weimar state FEs'),las=2,tick=F, cex.axis=.8)
points(coeffs,rev(pos),pch=20)
for (i in 1:4){
  segments(coeffs[i]-1.96*sts[i],rev(pos)[i],
           coeffs[i]+1.96*sts[i])
  segments(coeffs[i]-1.65*sts[i],rev(pos)[i],
           coeffs[i]+1.65*sts[i], lwd=2)
}
dev.off()

#### AFD + NPD 70km
coeffs <- c(-.171, cfs_currentstates_1st[4], cfs_weimar_centroid[4], cfs_weimar_contemp[4])
sts <- c(.066, ses_currentstates_1st[4], ses_weimar_centroid[4], ses_weimar_contemp[4])
pdf('elections_plot4.pdf', width=4.3, height=4.2)
par(mar=c(2.5,2,1,2))
plot(seq(-.35,.1, length.out = 2), c(.5,4.5), type='n',yaxt='n', ylab='',
     xlab='',cex.lab=.9)
abline(v=0,lty=2)
pos <- 1:4
points(coeffs,rev(pos),pch=20)
for (i in 1:4){
  segments(coeffs[i]-1.96*sts[i],rev(pos)[i],
           coeffs[i]+1.96*sts[i])
  segments(coeffs[i]-1.65*sts[i],rev(pos)[i],
           coeffs[i]+1.65*sts[i], lwd=2)
}
dev.off()



####
#### Approach 4: Weimar states based on areal overlap
####
# a) AFD, full sample
mod01_a <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateArea), elect_data)
mod01_b <- lm(I(AfDshare - coef(mod01_a)["foreign_rate"]*foreign_rate 
                - coef(mod01_a)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateArea), elect_data)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data[sample(1:nrow(elect_data), replace=TRUE),]
  boot1 <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateArea), databoot)
  boot2 <- lm(I(AfDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateArea), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs1 <- apply(gboot,2,sd)

# b) AFD, 70km
mod01_c <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateArea), elect_data70)
mod01_d <- lm(I(AfDshare - coef(mod01_c)["foreign_rate"]*foreign_rate 
                - coef(mod01_c)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateArea), elect_data70)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data70[sample(1:nrow(elect_data70), replace=TRUE),]
  boot1 <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateArea), databoot)
  boot2 <- lm(I(AfDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateArea), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs2 <- apply(gboot,2,sd)

# c) AFD+NPD, full sample
mod01_e <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateArea), elect_data)
mod01_f <- lm(I(AfDNPDshare - coef(mod01_e)["foreign_rate"]*foreign_rate 
                - coef(mod01_e)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateArea), elect_data)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data[sample(1:nrow(elect_data), replace=TRUE),]
  boot1 <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateArea), databoot)
  boot2 <- lm(I(AfDNPDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateArea), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs3 <- apply(gboot,2,sd)

# d) AFD+NPD, 70km
mod01_g <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateArea), elect_data70)
mod01_h <- lm(I(AfDNPDshare - coef(mod01_g)["foreign_rate"]*foreign_rate 
                - coef(mod01_g)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateArea), elect_data70)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data70[sample(1:nrow(elect_data70), replace=TRUE),]
  boot1 <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateArea), databoot)
  boot2 <- lm(I(AfDNPDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateArea), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs4 <- apply(gboot,2,sd)

###
### TABLE DA.2
###
stargazer(mod01_b, mod01_d, mod01_f, mod01_h,
          se=list(SEs1,SEs2,SEs3,SEs4),
          star.cutoffs = c(0.05, 0.01, NA), 
          column.sep.width="1pt", digits=3, no.space=T,
          dep.var.labels.include = FALSE)

cfs_weimar_area <- c(coef(mod01_b)[2], coef(mod01_d)[2], coef(mod01_f)[2], coef(mod01_h)[2])
ses_weimar_area <- c(SEs1[2],SEs2[2],SEs3[2],SEs4[2])



####
#### Approach 5: Weimar states based on centroids + non-matched Saarland areas as separate state
####
# a) AFD, full sample
mod01_a <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid2), elect_data)
mod01_b <- lm(I(AfDshare - coef(mod01_a)["foreign_rate"]*foreign_rate 
                - coef(mod01_a)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid2), elect_data)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data[sample(1:nrow(elect_data), replace=TRUE),]
  boot1 <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid2), databoot)
  boot2 <- lm(I(AfDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid2), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs1 <- apply(gboot,2,sd)

# b) AFD, 70km
mod01_c <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid2), elect_data70)
mod01_d <- lm(I(AfDshare - coef(mod01_c)["foreign_rate"]*foreign_rate 
                - coef(mod01_c)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid2), elect_data70)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data70[sample(1:nrow(elect_data70), replace=TRUE),]
  boot1 <- lm(AfDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid2), databoot)
  boot2 <- lm(I(AfDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid2), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs2 <- apply(gboot,2,sd)

# c) AFD+NPD, full sample
mod01_e <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid2), elect_data)
mod01_f <- lm(I(AfDNPDshare - coef(mod01_e)["foreign_rate"]*foreign_rate 
                - coef(mod01_e)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid2), elect_data)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data[sample(1:nrow(elect_data), replace=TRUE),]
  boot1 <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid2), databoot)
  boot2 <- lm(I(AfDNPDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid2), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs3 <- apply(gboot,2,sd)

# d) AFD+NPD, 70km
mod01_g <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid2), elect_data70)
mod01_h <- lm(I(AfDNPDshare - coef(mod01_g)["foreign_rate"]*foreign_rate 
                - coef(mod01_g)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid2), elect_data70)
boots <- 5000
set.seed(12435)
gboot <- matrix(NA, nrow=boots, ncol=6)
for(b in 1:boots){
  databoot <- elect_data70[sample(1:nrow(elect_data70), replace=TRUE),]
  boot1 <- lm(AfDNPDshare ~ distance2 + vshare33 + unemp33 + pop25 + prop_juden + unemp_rate + foreign_rate + as.factor(WeimarStateCentroid2), databoot)
  boot2 <- lm(I(AfDNPDshare - coef(boot1)["foreign_rate"]*foreign_rate - coef(boot1)["unemp_rate"]*unemp_rate) ~ distance2 + vshare33
              + unemp33 + pop25 + prop_juden + as.factor(WeimarStateCentroid2), databoot)
  gboot[b,] <- coef(boot2)[1:6]
  if (b%%250==0) print(b)
}
SEs4 <- apply(gboot,2,sd)

###
### TABLE DA.1
###
stargazer(mod01_b, mod01_d, mod01_f, mod01_h,
          se=list(SEs1,SEs2,SEs3,SEs4),
          star.cutoffs = c(0.05, 0.01, NA), 
          column.sep.width="1pt", digits=3, no.space=T,
          dep.var.labels.include = FALSE)

cfs_weimar_centroid2 <- c(coef(mod01_b)[2], coef(mod01_d)[2], coef(mod01_f)[2], coef(mod01_h)[2])
ses_weimar_centroid2 <- c(SEs1[2],SEs2[2],SEs3[2],SEs4[2])


####
#### Plots for FIGURE SM3.1
####

#### AFD full sample
coeffs <- c(cfs_weimar_centroid[1], cfs_weimar_centroid2[1], cfs_weimar_area[1])
sts <- c(ses_weimar_centroid[1], ses_weimar_centroid2[1], ses_weimar_area[1])
pdf('elections_weimaralt1.pdf', width=6.25, height=4.2)
par(mar=c(2.5,8,1,2))
plot(seq(-.15,.05, length.out = 2), c(.5,3.5), type='n',yaxt='n', ylab='',
     xlab='',cex.lab=.9)
abline(v=0,lty=2)
pos <- 1:3
axis(side= 2,at=rev(pos),labels=c('Centroid interpolation',
                                  'Centroid interpolation\n+ Saarland',
                                  'Area interpolation'),las=2,tick=F, cex.axis=.8)
points(coeffs,rev(pos),pch=20)
for (i in 1:3){
  segments(coeffs[i]-1.96*sts[i],rev(pos)[i],
           coeffs[i]+1.96*sts[i])
  segments(coeffs[i]-1.65*sts[i],rev(pos)[i],
           coeffs[i]+1.65*sts[i], lwd=2)
}
dev.off()

#### AFD 70km
coeffs <- c(cfs_weimar_centroid[2], cfs_weimar_centroid2[2], cfs_weimar_area[2])
sts <- c(ses_weimar_centroid[2], ses_weimar_centroid2[2], ses_weimar_area[2])
pdf('elections_weimaralt2.pdf', width=4.3, height=4.2)
par(mar=c(2.5,2,1,2))
plot(seq(-.35,.1, length.out = 2), c(.5,3.5), type='n',yaxt='n', ylab='',
     xlab='',cex.lab=.9)
abline(v=0,lty=2)
pos <- 1:3
points(coeffs,rev(pos),pch=20)
for (i in 1:3){
  segments(coeffs[i]-1.96*sts[i],rev(pos)[i],
           coeffs[i]+1.96*sts[i])
  segments(coeffs[i]-1.65*sts[i],rev(pos)[i],
           coeffs[i]+1.65*sts[i], lwd=2)
}
dev.off()

#### AFD + NPD full sample
coeffs <- c(cfs_weimar_centroid[3], cfs_weimar_centroid2[3], cfs_weimar_area[3])
sts <- c(ses_weimar_centroid[3], ses_weimar_centroid2[3], ses_weimar_area[3])
pdf('elections_weimaralt3.pdf', width=6.25, height=4.2)
par(mar=c(2.5,8,1,2))
plot(seq(-.15,.05, length.out = 2), c(.5,3.5), type='n',yaxt='n', ylab='',
     xlab='',cex.lab=.9)
abline(v=0,lty=2)
pos <- 1:3
axis(side= 2,at=rev(pos),labels=c('Centroid interpolation',
                                  'Centroid interpolation\n+ Saarland',
                                  'Area interpolation'),las=2,tick=F, cex.axis=.8)
points(coeffs,rev(pos),pch=20)
for (i in 1:3){
  segments(coeffs[i]-1.96*sts[i],rev(pos)[i],
           coeffs[i]+1.96*sts[i])
  segments(coeffs[i]-1.65*sts[i],rev(pos)[i],
           coeffs[i]+1.65*sts[i], lwd=2)
}
dev.off()

#### AFD + NPD 70km
coeffs <- c(cfs_weimar_centroid[4], cfs_weimar_centroid2[4], cfs_weimar_area[4])
sts <- c(ses_weimar_centroid[4], ses_weimar_centroid2[4], ses_weimar_area[4])
pdf('elections_weimaralt4.pdf', width=4.3, height=4.2)
par(mar=c(2.5,2,1,2))
plot(seq(-.35,.1, length.out = 2), c(.5,3.5), type='n',yaxt='n', ylab='',
     xlab='',cex.lab=.9)
abline(v=0,lty=2)
pos <- 1:3
points(coeffs,rev(pos),pch=20)
for (i in 1:3){
  segments(coeffs[i]-1.96*sts[i],rev(pos)[i],
           coeffs[i]+1.96*sts[i])
  segments(coeffs[i]-1.65*sts[i],rev(pos)[i],
           coeffs[i]+1.65*sts[i], lwd=2)
}
dev.off()
